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SUMMARY 

A semi -numerical method is reviewed for solving a set of coupled partial 
differential equations subject to mixed and possibly coupled boundary condi- 
tions. The line method of analysis is applied to the Navier-Cauchy equations 
of elastic and elasto-plastic equilibrium to calculate the displacement distri- 
butions in various, simple geometry bodies containing cracks. The application 
of this method to the appropriate field equations leads to coupled sets of 
simultaneous ordinary differential equations whose solutions are obtained along 
sets of lines in a discretized region. When decoupling of the equations and 
their boundary conditions is not possible, the use of a successive approxima- 
tion procedure permits the analytical solution of the resulting ordinary dif- 
ferential equations. 

The use of this method is illustrated by reviewing and presenting selected 
solutions of mixed boundary value problems in three-dimensional fracture 
mechanics. These solutions are of great importance in fracture toughness test- 
ing, where accurate stress and displacement distributions are required for the 
calculation of certain fracture parameters. Computations obtained for typical 
flawed specimens include that for elastic as well as ej as tor-plastic response. 
Problems in both Cartesian and cylindrical coordinate systems are included. 
Results are summarized for a finite geometry rectangular bar with a central 
through- the- thickness or rectangular surface crack under remote uniaxial ten- 
sion. In addition, stress and displacement distributions are reviewed for 
finite circular bars with embedded penny-shaped cracks, and rods with external 
annular or ring cracks under opening mode tension. 

The results obtained show that the method of lines presents a systematic 
approach to the solution of some three-dimensional mechanics problems with 
arbitrary boundary conditions. The advantage of this method over other numeri- 
cal solutions is that good results are obtained even from the use of a rela- 
tively coarse grid. 


INTRODUCTION 

The main goal of fracture mechanics is the prediction of the load at which 
a structure weakened by a crack will fail. Knowledge of the stress and dis- 
placement fields near the crack tip is of fundamental importance in evaluating 
this load at failure. Because of the geometric singularity associated with 
any crack type problem, there is almost no possibility of a simple closed form 
type of solution. For this reason, three-dimensional elastic solutions have 
been obtained only for a restricted class of problems. Furthermore, the calcu- 
lation of stress and strain distributions in elasto-plastic/work hardening 
materials containing inherent crack-like flaws is a nonlinear and three- 
dimensional problem. Due to the finite boundary effect and the nonlinearity 



of the material response, solutions in existence are obtained almost exclu- 
sively through numerical computer methods of continuum mechanics. Notable 
among these are the finite element method (refs. 1 and 2), the finite differ- 
ence method (ref. 3), and the boundary integral equation method (ref. 4). 

These methods are useful in solving either elastic or elasto-plastic fracture 
mechanics problems; it is known, however that practical problems usually 
require a very large amount of data storage and computation time. 

An alternate semi-analytical method suitable for the solution of crack 
problems is the line method of analysis. Successful application of this method 
to finite geometry solids containing cracks has been demonstrated recently for 
both elastic (ref. 5) and elasto-plastic (ref. 6) problems. Although the con- 
cept of the line method for solving partial differential equations is not new 
(refs. 7 and 8), its application in structural analysis has been limited to 
simple examples (ref. 9). By far the most common approach to fracture problems 
has been the finite element method, and it is the purpose of this paper to 
review a simple systematic, alternate method, the method of lines (MOL) for 
these problems. 

The line method lies midway between completely analytical and discretized 
numerical methods. The basis of this technique is the substitution of finite 
differences for the derivatives with respect to all the independent variables 
except one for which the derivatives are retained. This approach replaces a 
given partial differential equation with a system of ordinary differential 
equations whose solutions can then be obtained, at least in some cases, by ana- 
lytic methods. These equations describe the dependent variable along lines 
which are parallel to the coordinate in whose direction the derivatives were 
retained. Application of the line method is most useful when the resulting 
ordinary differential equations are linear and have constant coefficients 
(ref. 10). 

An inherent advantage of the line method over other numerical methods is 
that good results are obtained from the use of relatively coarse grids. This 
use of a coarse grid is permissible because parts of the solutions are obtained 
in terms of continuous functions. It is known that MOL methods tend to keep 
the advantages and discard the disadvantages of both the analytical and grid 
methods, thereby leading to accurate solutions with minimum computation times. 
Anot her a dvantage of this technique is that the displacements and the normal 
stress and strain fields are obtained with the same degree of accuracy at each 
node point. The disadvantage of MOL, on the other hand, is that it tends to 
become numerically unstable as the number of dividing lines increases and the 
finite difference strip size becomes too small (refs. 8, 11 and 12). To real- 
ize a very fine space discretization with t his method would require word length 
with much larger number of bits, leading to excessive requirements on computer 
resources. Current research emphasis in MOL solution methods is to overcome 
this problem in engineering applications (ref. 13). 


GOVERNING EQUATIONS AND MOL FORMULATION 

It is assumed, for simplicity of this presentation, that the material is 
homogeneous, isotropic and that the deformations are quasi-static and small. 

The structure is assumed to be elastic first and the elastic solution is taken 
to be known before the incipient loading is applied. As loading gradually 
increases, the structure becomes elasto-plastic and the governing equations are 
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written in terms of displacement increments. Using the standard summation con- 
vention, the Navier equations for the elastic problem in terms of displace- 
ments, u; are 


U i » j j 



i , j = 1,2,3 


( 1 ) 


and for the elasto-plastic regime, the displacement increments, duj , can be 
obtained from 


du. 


i.ll 



+ 3S. . 
iJ 



( 2 ) 


where the body forces are assumed to be zero, dep is the effective plastic 
strain increment, Sn is the stress deviator tensor and a e is the equivalent 
stress. In the plastic region the von Mises yield condition and the associated 
Prandl-Reuss flow rule is taken to prevail. The incremental stress-strain 
relations are obtained as (ref. 6), 
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(3) 


where v,G,E are the conventional elastic properties, ojj is the Kronecker 
delta and ajj are the stresses. 

In order to solve equations (1) or (2), we apply MOL and reduce these 
equations to systems of simultaneous, ordinary differential equations. For 
problems in Cartesian coordinates, equations (1) can be written as 

* (r^XD ■ « 

* feXg) - o <•> 

v2 » * (rH;)(S) - 0 (6) 


where the dilatation is 


3u 3v 3w 
e “ 9x + 3y + 3z 


( 7 ) 


and u, v, and w are the x-directional , y-directional and z-directional dis- 
placements, respectively. For a finite geometry solid with rectangular bounda- 
ries, we construct three sets of parallel lines (fig. 1(a)). Each set of lines 
is parallel to one of the coordinate axes and thus perpendicular to the corre- 
sponding coordinate plane. An approximate solution of equation (4) can then 
be obtained by developing solutions of ordinary differential equations along 
the x-directional lines. As seen in the figure, there are a total of 
il - NY x NZ such lines where NY is the number of lines along the y-direc- 
tion and NZ is the number of lines along the z-direction in a given plane, 
respectively. We define the displacements along these lines as ui,U 2 ,. . . , 
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(a) Three sets of lines parallel to x-, y - and z-coordinates 
and perpendicular to corresponding coordinate planes. 



(b) Set of interior lines parallel to x-coordinate. 

FIGURE 1. - SETS OF LINES PARALLEL TO CARTESIAN COORDINATES. 


uq,. The derivatives of the y-directional displacements on these lines with 
respect to y are defined as v'l^.v'^.. • -.v'Iq,, and the derivatives of 
the z-directional displacements with respect to z are defined as 

w* 1 1 ,w 1 1 2 w'ljj,. These displacements and derivatives can then be 

regarded as functions of x only since they are variables on x-directional 
lines. When these definitions are used the ordinary differential equation 
along a generic line ij (a double subscript is used here for simplicity of 
writing and the subscripts obviously are not related to those in the previous 
equations) in figure 1, using central differences with truncation errors of 
0(h 2 ) , may be written as 
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Similar differential equations are obtained along the other x-directional 
lines. The set of 2 second order differential equations represented by equa- 
tion (8) can be reduced to a set of 22 first order differential equations by 
treating derivatives of the u's as an additional set of 2 unknowns. The 
resulting equations can now be written as a single first order matrix differen- 
tial equation 

^ - A t U + R(x) ~ (10) 

where U and R are column matrices of 22 elements each and A* is 22 x 22 
matrix of coefficients. In a similar manner to solve the other two Navier 
equations for the elastic problem, we construct ordinary differential equations 
along the y- and z-directional lines, respectively. These equations are also 
expressed in an analogous form to equation (10); they are 

§ - v + S(y > (“I 

S " A 3 W + ( 12 ) 

Equations (10) to (12) are linear, first-order, ordinary differential equa- 
tions. They are, however, not independent, but are coupled through the vectors 
R, S and T. 

Noting that a second order ordinary differential equation can satisfy only 
a total of two boundary conditions and since three-dimensional elasticity prob- 
lems have three conditions at every point of the bounding surface, the shear 
stress boundary data must be incorporated into the differential equations of 
the surface lines. The applications of the specified shear conditions permits 
the use of a single layer of boundary image lines when surface line differen- 
tial equations are constructed. 
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For an elasto-plastic solid the governing differential equations for dis- 
placement increments and the incremental stress-displacement relations are 
found in reference 6. The x-directional displacement increments, in an analo- 
gous manner to equation (10), can be obtained from 

^(dU) =■ A x dU + dR(x) (13) 

where the coupling vector dR(x) contains mixed derivative terms for elastic 
and plastic regions in addition to terms involving the ratio of de p /a e . 

The system of ordinary differential equation (10) can be solved by any of 
a number of standard techniques. The method employed in references 5, 6, and 9 
is the Peano-Baker method of integration. The solution can be written as 

U(x) = e^UfO) + e^ e ^RfnJdn (14) 

J 0 

where U(0) is the initial value vector determined from the boundary conditions 

Ax 

and the matrizant e 1 is generally evaluated by its matrix series. For 
larger values of x, when convergence becomes slow, additive formulas may be 
used. In addition, similarity transformations can be used to diagonalize the 
coefficient matrix Aj. It should be noted that, in general, the matrix A^ 
is a function of Poisson's ratio and the coordinate finite difference incre- 
ments. Uniform line spacing in the three coordinate directions makes closed 
form diagonal izat ion of A^ possible. However, refinement of the mesh with 
uniform line spacing rapidly increases the required computer time and storage 
as well as raises the probability of numerical difficulties in the matrix expo- 
nential power series computations. Consequently, variable mesh spacing is 
recommended as one method of obtaining improved answers without an increase in 
problem size. 

Most of the previously obtained MOL solutions in fracture mechanics 
involved the use of finite difference formulas with truncation errors of 
0(h 2 ). Recent work by Mendelson and Alam (ref. 13), uses higher order finite 
difference approximations as an alternate method of obtaining more accurate 
results. These five point finite difference approximations for the first and 
second y-derivatives of a function f(x,y,z) at (x,y,z) can be written as 
(ref. 8): 
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x,y,z 


f (x,y + h y ,z - f (x,y - h y , z ) 


2h 


f (x,y + 2h y ,z) - f (x,y - 2h y , z ) 


4h 




w 

^3y 2 


x.y.z 


f(x,y + h ,z) + f(x,y - h ,z) - 2f(x,y,z) 
— — 


f (x,y + 2h , z) + f (x,y - 2h ,z) - 2f (x,y,z) 

J J - — 


4h: 


°W) < i5) 


Since the evaluation of the matrix exponential power series becomes increas- 
ingly difficult with the coefficients obtained through the use of equa- 
tions (15), a recurrence relations method is used to solve the resulting sys- 
tems of ordinary differential equations. An error analysis in reference 8 
indicates that approximately six times as many dividing lines must be used 
with 0(h 2 ) approximations to get equivalent accuracy to that obtained when 
equations (15) are used. 


The solution of equation (10) by recurrence relations can be obtained by 
taking N equal intervals along the x-axis, each having a length of h ra . 

Then by using finite differences for (dU/dx) m and average values of 
(A^U + R) m , the following linear recurrence formula will' be obtained 
(ref. 13), expressing U m in terms of U m _i: 


U 


m 


LU , + M 
m m-1 m 


(16) 


where Ln, is a known function of h m and A^ while M m will depend on R 
in addition to h m and Aj. We can also express U m in terms of Uj by 
repeated application of equation (16), leading to 


U = D IL + 
m ml 


m 


(17) 


where Djj and F m are known functions of L m and M m . By suitable parti- 
tioning of D and F matrices at the boundaries, we can use given boundary 
data at the last station, U n , to calculate unknown elements of Uj, where 
n = N + l. The advantage of using equation (17) to calculate U m , 
m = 1,2,. . .,n, is that the coefficient matrix A^ has no limitation on its 
format or on the arrangement of its elements. The interval h m can be 
decreased to any fraction of h x , the initially established finite difference 
increment obtained from the application of MOL. Results of test problems indi- 
cate that two or three subintervals are adequate for the solution of a typical 
problem. 

In a similar manner, for problems in cylindrical coordinates, 
equations (1) can be written as 
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(3e/3r) + (1 - 2v)[(r - r ^)u - 2r *(3v/30)] - 0 


(18 


r 1 ( 3e/3©) + (1 - 2v)[(V 2 - r 2 )v + 2r 2 (3u/30)] = 0 
(3e/3z) + (1 - 2 v)V 2 w = 0 

where the dilatation is 

e - (3u/3r) + r 1 (3v/30) + (u/r) + (3w/3z). 


(19 

(20 



FIGURE 2. - SETS OF LINES IN DIRECTION OF CYLINDRICAL 
COORDINATES. 


For a finite geometry body with circular boundaries, we construct three 
sets of parallel lines in the direction of the coordinates as shown in fig- 
ure 2. Approximate solutions of the above equations can then be obtained by 

developing solutions of ordinary differential equations along the radial, cir- 
cumferential, and axial lines, respectively. For equation (18), we define the 
displacements along the radial lines as ui,U 2 ,. . . .uj,. The derivatives of 

the circumferential displacements on these lines with respect to 0 are 

defined as v ’ j ^ , v ' 1 2 v 'lsi> the derivatives of the axial displace- 

ments with respect to z are defined as w'l^.w'^,. • . ,w'|g,. These dis- 
placements and derivatives can then be regarded as functions of the radius 
only, since they are variables on radial lines. If these definitions are used, 
the ordinary differential equation along a generic, radial line ij of fig- 
ure 2 may be written as 
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l d u ii 1 du ii 
dr 2 + r dr 


u ij 


(1 - 2y ) 
+ 2(1 - v) 


O Z W 


+ -r(u. . . + u. . , ) 

h 2 l i,j+l i.J-l' 

z 


f. . (r) 

* zir^r ' 0 (21) 


where 


f i ; ( r ) - [(4v - 3)/r 2 ]v'|.. + r 1 (dv’/dr)| i . + (dw'/dr)|.. 


( 22 ) 


and 


v' = dv/d© 
w' - dw/dz 

Similar differential equations are obtained along the other radial lines. 

Since each equation has the terms of the displacements on the surrounding 
lines, these equations constitute a system of ordinary differential equations 
for the displacements ui,U 2 ,. . . ,u&. 

Noting that a second-order differential equation can satisfy only a total 
of two boundary conditions, the shear stress boundary data must again be incor- 
porated into the surface line differential equations. For the first radial 
line of figure 2, the use of zero shear stress boundary conditions in the 
radial direction on the r,z and r,0 coordinate planes gives, respectively, 
the following imaginary radial line displacements: 

u i©e - u 2 + 2 V (dv/dr) l 1 - 2h e v| i 

u 16n - U N6+1 + 2h z ( dw/dr, ll (23) 


Equations (23) must then be used in the application of central difference 
approximations when the ordinary differential equation for the first radial 
line is generated. Additional details on the construction of these equations 
can be found in reference 13. 

As before the system of & second order equations given by equation (21) 
can be written as a single first order matrix equation 

dU/dr - A r (r)U + R(r) (24) 


where in this case the elements of U are defined by 


"i ■ ru i ■ 


U 2 - rU 2 , 


U l ' rU H 


u 


8.+1 


r 1 [d(ru 0 )/dr] , U a+2 = r 1 [d(ru 9 )/dr] , . . .,, 


(25) 

U 2!l - r -1 [d(ru a )/dr] 
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and the elements of the coefficient matrix A r (r) are no longer constant, but 
are functions of r. The corresponding equations for V and ff are the same 
as equations (11) and (12) with y replaced by 0, the coefficient matrices 
being constant. 

The systems of ordinary differential equations for the circumferential 
and axial displacements are solved in the same manner as equations (11) and 
(12), that is, by using an analogous solution to equation (14). For the case 
of equation (24) the solution is more complicated, since we no longer have con- 
stant coefficients. The solution may be expressed as (ref. 5). 


U(r) = G(A r )U(r 0 ) + Q(A f ) 


(Q(A r )] 1 R(n)dn 


(26) 


where rg is the initial value of r which may or may not be zero, 
matrizant Q(A r ) is given by the infinite matrix integral series 


The 


r-r 


Q(A r ) - I + 


nr 


A f (n 1 )dn 1 + 


n 2 


V n 2 )dn 2 


A r (n 1 )dn 1 


r r 

r 3 r 

A r (n 3 ) d n 3 + 

A r (n 2 )dn 2 

J r n 

J 


A^(n 1 )dn 1 (27) 


which of course becomes very difficult to evaluate. However, by substituting 
the matrix A r into equation (27), it can be seen by inspection, (ref. 5), 
that Q(A r ) can be partitioned into four submatrices which satisfy the follow- 
ing equations: 


wi th 


r 1 (dQ 11 /dr) = Q 21 
r(dfl 21 /dr) - Q n K r 

r _ 1 (dn i 2 /dr) = Q 22 
r(dfl 22 /d r) = Q 12 K r 


n il'r=r 0 " Q 22 lr=r 0 " 1 


^12'r=rg = ^21'r=rg = 0 


( 28 ) 
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and K r is the submatrix of A r obtained from the identity: 

0 rl 
r _1 K r 0 

Instead of evaluating the matrizant Q(A r ) by means of equation (27), it 
can be evaluated by solving equations (28). This was done for the examples 
presented herein, using a single step Runge-Kutta integration method for solv- 
ing equations (28). 

Since equations (10) to (12) and their boundary conditions are highly cou- 
pled, it is generally impossible to directly evaluate their solutions. Thus, a 
successive approximation procedure must be employed where assumed values must 
be used initially for the required unknowns. Once the displacement field in 
the body has been calculated and the successive approximation procedure has 
converged, the normal stress distributions can be obtained directly by using 
the stress-displacement relations. The shear stresses, however, can be evalu- 
ated only through finite difference approximations for the required displace- 
ment gradients. 

All of the MOL work in three-dimensional fracture mechanics has been done 
using double precision arithmetic. With larger word sizes, 128 bits or 
greater, improved results can be obtained or more lines can be used. Typical 
problems on third generation computers can usually be handled with 100 to 200 
lines in each direction using Cartesian coordinates, and up to 20 lines in each 
direction using cylindrical coordinates. Corresponding CPU times are of the 
order of 3 min for the elastic solution and 25 min for the entire elasto- 
plastic problem (ref. 6). Most of the computer time is spent in decoupling the 
three systems of simultaneous ordinary differential equations which arise in a 
general problem. Decoupling is done by a successive approximation procedure. 
With cyclic resubstitution of the obtained solutions into the coupling vectors 
and the boundary conditions good numerical convergence behavior was observed. 

Elastic solutions have been obtained for typical fracture test specimen 
geometries such as the central crack, single edge crack, double edge crack and 
rectangular surface crack problems, all with uniform tensile loadings normal 
to the crack plane. In cylindrical coordinates, the problem of an embedded 
penny-shaped crack and the externally cracked circular cylinder in tension 
have been treated. Although shear and torsional ly loaded specimens have not 
been analyzed previously by MOL, the necessary boundary conditions for these 
problems can be imposed without difficulties. Recently, the thumbnail crack 
problem was investigated in connection with applying MOL to curved crack bound- 
aries (ref. 13). In addition, it seems that the common three-point bend speci- 
men could also be analyzed in a systematic manner using this method. 

Elasto-plastic solutions of certain crack problems have also been 
obtained using the incremental displacement formulation in connection with MOL 
(ref. 6). The nonlinear response of finite length cylinders with external 
annular cracks and a finite thickness rectangular plate containing a through- 
thickness central crack under uniaxial tension were studied in detail. In 
addition to the stresses and displacements, fracture mechanics parameters such 
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as the stress intensity factor, the J-integrai and the load versus load point 
displacement plots were determined. Finally, the thumbnail crack problem was 
also studied, permitting small scale nonlinear material response. 


STRESS INTENSITY AND STRAIN ENERGY DENSITY FACTORS 

Since the application of boundary conditions for MOL allows the crack tip 
to remain between two successive node points, the exact location of crack tip 
together with the determination of K values for the elastic case is done by 
the first two terms in the Williams eigenfunction expansion. The crack face 
displacement and maximum normal stress near the crack front are used to find 
the coefficients in the two-term expansion. Assuming Cartesian coordinates and 
that y = 0 is the crack plane and that the crack is under normal tensile load- 
ing, we have 


v = otK, 


yEp ♦ ^VUTTF) 3 


o = K r 

y i 


i— .?VS-T7 


LV2ir(R - r) I 


(29) 


(30) 


where a is a function of Poisson's ratio, the stress singularity is assumed 
to be -1/2, r is the crack edge position correction, v is the crack opening 
displacement, oy is the maximum normal stress and Kj and are the 

mode I Williams expansion coefficients. Using displacement data from three 
adjacent nodes to the crack edge in equation (29), values of aKi , L j /Kj and 
r are calculated for each value of z, with R also measured from the half- 
way point between nodes specifying boundary stresses and displacements, respec- 
tively. Substituting values of L j /K j and r into equation (30), we can 
calculate Kj as a function of the corrected crack edge distance, p » R - r. 
Note that a would be equal to 3.56 for the plane strain case and 4.0 for the 
plane stress case with v = 1/3. Another approach to calculate Kj is to 
first determine the Ji integral and then use the linear elastic J[ - Kj 
relation of the form (ref. 6) 



Linear fracture mechanics technology assumes then that the crack will 
propagate if Kj reaches its critical value Kj c , usually called its fracture 
toughness. It should be noted that the K-concept is restricted to symmetric 
systems with the applied loads perpendicular to the crack plane and the crack 
propagating in a self-similar manner. In general, a complete description of 
the crack border stress field requires three stress intensity factors, and a 
mixed mode fracture criterion is needed to predict failure. To this end Sih 
(ref. 14) has defined a strain energy density factor, S, as 

S - + ^ a i2^1^2 + a 22^2 + a 33^3 
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where the coefficients j ( i , j - 1,2) depend on the material constants and 
the angles 6 and u. Consistent with equation (32), the local stresses near 
the crack tip are of the form (ref. 14), 


i e . e . 39A 

<, - -== ™ 2V 1 - Sln 2~ sln 2 ~) 

V2pcosco 


- 2 sin f-f2 - cos f-cos + 0(1) (33) 

V2pCOS(d 

Note that ki = Kj/\/ir and the coefficients a^j are then given by equations 
of the form 

16G cos to = (3 - 4v - cos 9)(1 + cos ©) (34) 

As can be seen from equation (32), the stress intensity factors kj still play 
an important role in the fracture process. Hence, the correct determination 
of these factors is a necessary step in the safe design of any structure. In 
the strain energy density failure criterion, it is assumed that the minimum 
value of S yields the direction of crack initiation and that the critical 
value of S m j n , S c , determines incipient fracture and is an intrinsic material 
property independent of the loading conditions and crack configurations. 

Other multimode failure criteria have been proposed previously in the lit- 
erature and a brief description of each fracture theory along with the limited 
experimental data can be found in (refs. 15 and 16). In general, little dif- 
ference exists between theories predicting mode I, mode II interaction and com- 
bined damage crack propagation direction. 


NUMERICAL RESULTS 

A great amount of experimental work has been done in fracture mechanics 
using cracked specimens. Selected results for some common specimen geometries 
are shown in figures 3 to 16. Figure 3 shows a rectangular bar under normal 
tensile loading containing a traction free, through- thi ckness , central crack. 
The crack opening displacement at the middle and at the surface of the bar is 
plotted in figure 4, along with Raju's (ref. 17) finite element results. 

Stress intensity factor variations as a function of bar thickness are shown in 
figure 5. Variation of the constraint parameter 3, defined as the quantity 
<7 z /[v(a x + <jy) ] , along the plate thickness is shown in figure 6. Note that 
for plane strain 3-1 and for plane stress case it vanishes. Figure 7 shows 
a bar with uniform tension containing a rectangular surface crack (ref. 18). 
Surface crack opening displacement as a function of crack depth is shown in 
figure 8 while the variation of the maximum normal stress ay is shown in fig- 
ure 9 for a selected crack geometry. For the same rectangular surface crack 
problem, a plot of Kj along the crack periphery is shown in figure 10. The 
discretization of an externally cracked cylindrical fracture specimen is shown 
in figure 11. Crack face displacements for various crack lengths are plotted 
in figure 12 while the variation of Kj with crack length is shown in fig- 
ure 13. Figure 14 is analogous to figure 11 but with an embedded, central, 
penny-shaped crack. The dimensionless crack opening displacements for various 
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<al RECTANGULAR BAR WITH THROUGH -THICKNESS CENTRAL 
CRACK. 
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(b/ DISCRETIZED REGION OF RECTANGULAR BAR WITH 
THROUGH-THICKNESS CENTRAL CRACK. 

FIGURE 3. - RECTANGULAR BAR WITH THROUGH-THICKNESS CENTRAL 
CRACK UNDER UNIFORM TENSION. 
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FIGURE 5. - STRESS- INTENSITY FACTOR VARIATION AS A 
FUNCTION OF BAR THICKNESS FOR A CENTER-CRACKED REC- 
TANGULAR BAR UNDER UNIFORM TENSION. 
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FIGURE 6. - VARIATION OF 0 = O/lVtO* ♦ Oy) ) ON CRACK PLANE < y * 0. 
ACROSS PLATE THICKNESS. (c/W) « 0.517M. FOR IDEAL PLANE STRAIN 
CASE: P = 1, AND FOR PLANE SIRESS: 0 = 0. 



FIGURE 7. - BAR WITH RECTANGULAR SURFACE CRACK UNDER UNIFORM TENSION. 





FIGURE 8. - SURFACE CRACK OPENING DIS- 
PLACEMENT VARIATION AS A FUNCTION OF 
CRACK DEPTH FOR A RECTANGULAR BAR 
UNDER UNIFORM TENSION. 



(a> DIMENSIONLESS y-DIRECTIONAL NORMAL STRESS 
VARIATION ALONG BAR WIDTH. 


h\ 



(b) DIMENSIONLESS y-DIRECTIONAL NORMAL STRES5 
VARIATION ACROSS BAR THICKteSS. 

FIGURE 9. - DIMENSIONLESS y-DIRECTIONAL NORMAL 
STRESS DISTRIBUTION IN THE CRACK PLANE FOR A 
BAR UNDER UNIFORM TENSION CONTAINING A RECTAN- 
GULAR SURFACE CRACK, 
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FIGURE 15. - DIMENSIONLESS CRACK OPENING DISPLACE- 
MENTS FOR SOLID CYLINDRICAL BARS WITH PENNY-SHAPED 
CRACKS OF VARIOUS LENGTHS AND RADII. 



FIGURE IN. - SOLID CYLINDRICAL BAR WITH PENNY-SHAPED 
CRACK. 


FIGURE 16. - DIMENSIONLESS MAXIMUM CRACK OPENING AS 
FUNCTION OF CRACK TO CYLINDER RADIUS RATIO. 
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geometry cylinders are plotted in figures 15 and 16 and are compared to 
Sneddon's (ref. 19) closed form solution. It is obvious from these results 
that a variety of plots familiar to the fracture mechanics community can be 
constructed, since MOL methods give complete field solutions. 


CONCLUDING, REMARKS 

The line method of analysis is a practical approach for the solution of 
three-dimensional crack problems, at least for bodies with reasonably regular 
boundaries. Both elastic and inelastic solutions can be obtained. Just how 
efficient the method is or can be made is not fully established. It is known, 
however, that good results are obtained from the use of relatively coarse 
grids. Interestingly, displacements and normal stresses are determined with 
equal accuracy since numerical differentiation of the displacements is not 
required. Applications to curved boundaries, bending or shear modes of loading 
and variable mesh spacing are some of the current areas that need additional 
investigations. Furthermore, improvements in computer hardware, such as the 
use of super computers, could help overcome some of the presently experienced 
stability problems. 
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